---
title: "Replication and Extension_Gov2020"
author: "David Aboge, Nidhi Patel, Nic Wicaksono"
paper: "Trauma /& Turnout: The Political Consequences of Traumatic Events"
date: "Nov 2024
output:
---

# Setting up the Data {.tabset}

## Loading packages and data

```{r setup, message = FALSE, warning = FALSE, results ='hide'}

##### Packages #####

library(arm)
library(plyr)
library(boot)
library(descr)
library(ggplot2)
library(psych)
library(GPArotation)
library(dplyr)
library(RPMG)
library(foreign)
library(ggpubr)
library(ggeffects)
library(data.table)
library(lubridate)
library(haven)
library(panelView)
library(plm)
library(multiwayvcov) 
library(lmtest)
library(wfe)
library(PanelMatch)
library(DataCombine)
library(mi)
library(stringr)
library(lfe)
library(clubSandwich)
library(reghelper)
library(stargazer)
library(xtable)

setwd("I:/.shortcut-targets-by-id/1mzU-tavrkNyohejb8pBpzc1C6tOQCgbR/Replication Project/Marsh 2022 - Voting and Trauma/replication files") # set working directory
knitr::opts_knit$set(root.dir = 'I:/.shortcut-targets-by-id/1mzU-tavrkNyohejb8pBpzc1C6tOQCgbR/Replication Project/Marsh 2022 - Voting and Trauma/replication files')

```

# Load data

```{r}

load('cps.Rda') # individual-level CPS data
load('cps_lag.Rda') # individual-level CPS data w/ one lag
load('cps_2lag.Rda') # individual-level CPS data w/ two lags
load('discop.Rda') # county-level data
load('stanrob.Rda') # Stanford Mass Shootings data for robustness check

# refer to "cps.Rda" as cps3 in code
# refer to "cps_lag.Rda" as cps_lag in code
# refer to "cps_2lag.Rda" as cps_2lag in code
# refer to "discop.Rda" as discop in code
# refer to "stanrob.Rda" as stan2 in code

```


## Creating Incumbent Vote Share Variable

```{r, include = F} 

table(discop$inc_part, discop$year)
# -1 Republican, 1 Democrat, 0 if there is no incumbent president

discop$rul_part <- discop$inc_part
discop$rul_part[discop$year == 1988 | discop$year == 2008] <- -1
discop$rul_part[discop$year == 2000 | discop$year == 2016] <- 1
table(discop$rul_part, discop$year)
# recoding years with no incumbent

discop$rul_part[discop$rul_part == -1] <- 0 
# making Republican incumbents "0" for ease of analysis 

discop$rul_vote_share[discop$rul_part == 1] <- discop$votes_dem/discop$votes_tot
discop$rul_vote_share[discop$rul_part == 0] <- discop$votes_gop/discop$votes_tot
# creating incumbent vote share

```

## Note that some tables were further customized prior to be placed in the working paper. These were only about collapsing the various control variable coefficients into a single "Controls" or "Fataliites and Injuries" category, for presentational purposes. The code below will therefore generate more detailed results. 

## Generalized Two Way FEs Analysis 

```{r, include = F}

library(clubSandwich)

# arson

art_plm1 <- plm(rul_vote_share ~ tr_sh + perc_bl + pop_tot + income, # with controls
           data = discop,
           index = c('fips', 'year'),
           model = 'within',
           effect = 'twoways')

# shooting

sht_plm1 <- plm(rul_vote_share ~ tr_ar + perc_bl + pop_tot + income +
                  fatalities + injured, # with controls
           data = discop,
           index = c('fips', 'year'),
           model = 'within',
           effect = 'twoways')

# natural disasters

dist_plm2 <- plm(rul_vote_share ~ tr_dis2 + perc_bl + pop_tot + income, # with controls
           data = discop,
           index = c('fips', 'year'),
           model = 'within',
           effect = 'twoways')


# Table 1 for Paper 

stargazer(art_plm1, sht_plm1, dist_plm2, 
          title = "Effect of Traumatic Events on Ruling Party Vote Share, County Level", 
          dep.var.caption = "Incumbent Presidential Party Vote Share", 
          dep.var.labels = "", 
          covariate.labels = c("Black Arson", "Mass Shootings", "Natural Disasters", "Pct. Population, Black", "Total Population", "Median HH Income", "Fatalities, Shootings", "Injuries, Shootings"), 
          add.lines = list(c("County \\& Year Fixed Effects","$\\checkmark$", "$\\checkmark$", "$\\checkmark$"))
          )

```


## Adding Current Incumbent Party as a Control 

``` {r, eval = F} 


## With controls

# arson

ar_lag2 <- plm(rul_vote_share ~ rul_part + tr_ar +
                perc_bl + pop_tot + income,
              data = discop, 
              index = 'fips', 
              model = 'within', 
              effect = 'individual')

summary(ar_lag2) 

# shooting

sh_lag2 <- plm(rul_vote_share ~ rul_part + tr_sh +
                perc_bl + pop_tot + income + fatalities + injured,
              data = discop, 
              index = 'fips', 
              model = 'within', 
              effect = 'individual')

summary(sh_lag2) 

# natural disaster

nat_lag2 <- plm(rul_vote_share ~ rul_part + tr_dis2 +
                perc_bl + pop_tot + income,
              data = discop, 
              index = 'fips', 
              model = 'within', 
              effect = 'individual')

summary(nat_lag2) 

# Table X For Main Paper - ultimately scrapped

stargazer(ar_lag2, sh_lag2, nat_lag2, 
          title = "Effect of Traumatic Events on Presidential Ruling Party Vote Share, County Level Votes, with Incumbent Party Control Variable", 
          dep.var.caption = "Incumbent Presidential Party Vote Share", 
          dep.var.labels = "", 
          covariate.labels = c("Dem. Incumbent", "Black Arson",  "Mass Shootings", "Natural Disasters", "Pct. Population, Black", "Total Population", "Median HH Income", "Fatalities, Shootings", "Injuries, Shootings"),
          add.lines = list(c("County Fixed Effects", "$\\checkmark$", "$\\checkmark$", "$\\checkmark$"))
          )

```

## Adding Interaction Effects - Party and Traumatic Event

``` {r, include = F} 


## With controls

# arson - no observations for interaction effect due to data limitations

 ar_int <- plm(rul_vote_share ~ rul_part + tr_ar +
               perc_bl + pop_tot + income,
             data = discop, 
             index = 'fips', 
             model = 'within', 
              effect = 'individual')

summary(ar_int) 

# shooting

sh_int <- plm(rul_vote_share ~ rul_part*tr_sh +
                perc_bl + pop_tot + income + fatalities + injured,
              data = discop, 
              index = 'fips', 
              model = 'within', 
              effect = 'individual')

summary(sh_int)

# natural disaster

nat_int <- plm(rul_vote_share ~ rul_part*tr_dis2 +
                perc_bl + pop_tot + income,
              data = discop, 
              index = 'fips', 
              model = 'within', 
              effect = 'individual')

summary(nat_int) 

# Table 3 for the paper


stargazer(ar_int, sh_int, nat_int, 
          title = "Effect of Traumatic Events on Ruling Party Vote Share, County Level, with Incumbent Party Control Variable", 
          dep.var.caption = "Incumbent Presidential Party Vote Share", 
          dep.var.labels = "", 
          covariate.labels = c("Dem. Incumbent", "Black Arson",  "Mass Shootings", "Natural Disasters", "Pct. Population, Black", "Total Population", "Median HH Income", "Fatalities, Shootings", "Injuries, Shootings", "Dem. Incumbent * Shootings", "Dem. Incumbent * Disaster"),
          add.lines = list(c("County Fixed Effects", "$\\checkmark", "$\\checkmark", "$\\checkmark"))
         )


```

## Annex 1 - Check Using Cumulative Incidents Leading Up to Election Year

``` {r, include = F} 


## With controls


ar_int <- plm(rul_vote_share ~ rul_part*cum_ar +
                perc_bl + pop_tot + income,
              data = discop, 
              index = 'fips', 
              model = 'within', 
              effect = 'individual')

summary(ar_int) 

# shooting
# No difference, since it's one incident per county-year anyways 

sh_int <- plm(rul_vote_share ~ rul_part*cum_sh +
                perc_bl + pop_tot + income + fatalities + injured,
              data = discop, 
              index = 'fips', 
              model = 'within', 
              effect = 'individual')

summary(sh_int) 

# natural disaster

nat_int <- plm(rul_vote_share ~ rul_part*cum_dis +
                perc_bl + pop_tot + income,
              data = discop, 
              index = 'fips', 
              model = 'within', 
              effect = 'individual')

summary(nat_int) 

# Table 4


stargazer(ar_int, sh_int, nat_int, 
          title = "Effect of Traumatic Events on Ruling Party Vote Share, County Level, with Incumbent Party Control Variable", 
          dep.var.caption = "Incumbent Presidential Party Vote Share", 
          dep.var.labels = "", 
          covariate.labels = c("Incumbent Party,\\\\ 1 = Democratic", "Black Arson",  "Mass Shootings", "Natural Disasters", "Pct. Population, Black", "Total Population", "Median HH Income", "Fatalities, Shootings", "Injuries, Shootings", "Inc. Party*Shootings", "Inc. Party*Disaster"),
          add.lines = list(c("County Fixed Effects", "$\\checkmark", "$\\checkmark", "$\\checkmark"))
         )


```

## Annex 2 - Using Single Measure of Whether Some Incident Happened 

```{r} 

# Generating Indicator 

discop$tr_any <- 0 
discop$tr_any[discop$tr_sh == 1 | discop$tr_dis2 == 1 | discop$tr_ar == 1] <- 1


## Running the Model 

any_tr <- plm(rul_vote_share ~ rul_part*tr_any + perc_bl + pop_tot + income, 
              data = discop, 
              index = 'fips', 
              model = 'within',
              effect = 'individual')

summary(any_tr) 

# Table 5

stargazer(any_tr, 
          title = "Predicting Incumbent Party Vote Share, using Pooled Trauma Event Measure", 
          dep.var.caption = "Incumbent Presidential Party Vote Share", 
          dep.var.labels = "", 
          covariate.labels = c("Dem Incumbent", "Any Trauma Event", "Pct. Population, Black", "Total Population", "Median HH Income", "Dem. Incumbent * Disaster"),
          add.lines = list(c("County Fixed Effects", "$\\checkmark", "$\\checkmark", "$\\checkmark"))
         )




```


## Descriptives

``` {r, include = F} 

summary_table <- discop %>%
  filter(year %in% c(1980, 1984, 1988, 1992, 1996, 2000, 2004, 2008, 2012, 2016)) %>% # Filter for specific years
  group_by(year) %>%
  summarise(
    mean_rul_vote = mean(rul_vote_share, na.rm = TRUE), 
    sd_rul_vote = sd(rul_vote_share, na.rm = TRUE),
    mean_tr_ar = mean(tr_ar, na.rm = TRUE),
    sd_tr_ar = sd(tr_ar, na.rm = TRUE),
    mean_tr_sh = mean(tr_sh, na.rm = TRUE),
    sd_tr_sh = sd(tr_sh, na.rm = TRUE),
    mean_tr_dis2 = mean(tr_dis2, na.rm = TRUE),
    sd_tr_dis2 = sd(tr_dis2, na.rm = TRUE),
    mean_perc_bl = mean(perc_bl, na.rm = TRUE),
    sd_perc_bl = sd(perc_bl, na.rm = TRUE),
    mean_pop_tot = mean(pop_tot, na.rm = TRUE),
    sd_pop_tot = sd(pop_tot, na.rm = TRUE),
    mean_income = mean(income, na.rm = TRUE),
    sd_income = sd(income, na.rm = TRUE),
    mean_fatalities = mean(fatalities, na.rm = TRUE),
    sd_fatalities = sd(fatalities, na.rm = TRUE),
    mean_injured = mean(injured, na.rm = TRUE),
    sd_injured = sd(injured, na.rm = TRUE)
  )
print(summary_table)
latex_table <- xtable(summary_table)

print(latex_table, file = "summary_table.tex")

```


